Statistics of radiation at Josephson parametric resonance 
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Motivated by recent experiments, we study theoretically the full counting statistics of radiation 
emitted below the threshold of parametric resonance in a Josephson junction circuit. In contrast to 
most optical systems, a significant part of emitted radiation can be collected and converted to an 
output signal. This permits studying the correlations of the radiation. 

To quantify the correlations, we derive a closed expression for full counting statistics in the limit 
of long measurement times. We demonstrate that the statistics can be interpreted in terms of 
uncorrelated bursts each encompassing 2N photons, this accounts for the bunching of the photon 
pairs produced in course of the parametric resonance. We present the details of the burst rates. In 
addition, we study the time correlations within the bursts and discuss experimental signatures of 
the statistics deriving the frequency- resolved cross-correlations. 

PACS numbers: 74.50-|-r, 73.23Hk, 85.25Cp 



I. INTRODUCTION 

Parametric resonance^ is one of the most fundamen- 
tal and frequently applied non-linear phenomena. If a 
non-linear oscillator with the resonant frequency fio is 
a.c. driven at frequency 257 2f2o, a coherent reso- 
nant response at frequency Vt emerges provided the driv- 
ing amplitude exceeds an instability threshold set by the 
non-linear parameters of the oscillator. While the coher- 
ent response is absent below the threshold, the paramet- 
ric resonance is manifested there by enhanced fluctua- 
tions with frequencies close to 51. In the quantum realm 
[hVt ^ ksT with T the temperature), these fluctuations 
can be regarded as an emission of radiation. An elemen- 
tary radiation event is an emission of a pair of photons of 
the frequency ~ 51 caused by absorption of a single pho- 
ton of the frequency 2J1. In quantum optics, the corre- 
sponding phenomenon is called down-conversion^ since a 
single photon is converted into two. The down-conversion 
is a base of optical quantum-information applications"^. 
The phenomenon has been employed to produce squeezed 
states of lighti and pairs of quantum-entangled photons^. 

It seems natural to assume that the statistics of the 
radiation is that of uncorrelated elementary events, each 
event being the emission of a correlated/entangled pair. 
In most optical experiments, this assumption is correct 
and practical. However, it relies on the fact that only a 
minor fraction of emitted pairs is actually detected. De- 
tected events are separated by large time intervals and 
thus do not show any correlation. Recently, a set of pio- 
neering experimentsS has advanced quantum non-linear 
optics into the microwave frequency range. Thereby, 
atoms are replaced by superconducting qubits made us- 
ing Josephson junctions, and the radiation is confined to 
transmission lines and electrical oscillators. The latter 
represents a large technical advantage in comparison with 
an optical experiment due to the fact the radiation is not 
lost and concentrated. This enhances the non-linearities 
of the system. 

Very recently, accurate measurements of the radiation 



emitted by a dc voltage-biased Josephson junction em- 
bedded in a microwave resonator have been reported 
The Josephson generation frequency loj = 2eV/h can 
be tuned to double the resonant frequency, fulfilling the 
conditions of parametric resonance. Importantly, up to 
50% of the emitted radiation can be detected and the 
fiuctuations of the detector signal can be quantified as 
welli^ This motivated us to study the statistics of the 
radiation in this setup. The full photon counting statis- 
tics of the degenerate optical parametric oscillator has 
been addressed ir*S for a specific case when the driving 
frequency is precisely 2VIq. Since the observation of non- 
Poissonian features of these statistics requires collection 
efficiency not achievable in optical setups, this work has 
not attracted the attention it deserves. Let us note that 
the measurements of statistics do not require the detec- 
tor to be an actual counter giving the output signal in 
terms of discrete numbers of counts. A continuous de- 
tector output would suffice to quantify cumulants of the 
radiation intensity fluctuations and thereby characterize 
the statistics. 

In this paper, we revisit the full counting statistics 
(FCS) of radiation below the instability threshold bring- 
ing this to the context of Josephson circuit. This regime 
is interesting since despite the fact that the fleld corre- 
lations are entirely Gaussian under these conditions, the 
statistics are highly nontrivial. We restrict our attention 
to FCS in the limit of the long measurement times. We 
recover the results o& in a different conceptual framework 
that is directly based on the Keldysh-action treatment of 
dissipative Josephson dynamics. We extend the results 
to the case of an arbitrary mismatch between driving 
frequency and 2flo. We provide an interpretation of the 
statistics. In this interpretation, an elementary event is 
a radiation burst that encompasses correlated emission 
of N pairs, rather than an emission of a single pair. This 
is a manifestation of photon bunching. We outline the 
similarities with the results ^"^concerning the bunching in 
a single-photon regime. The rate of iV-burst does not di- 
verge upon approaching the threshold. However, larger 
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TV are favored closer to the threshold. This results in a di- 
vergence of the average radiation intensity and its higher 
moments. We support this interpretation by investigat- 
ing time correlations of the emission events. Further, we 
quantify the frequency-resolved fluctuations of the radi- 
ation. The correlations of the spectral-resolved intensity 
permit a relatively easy experimental observation and we 
present several relevant formulas to facilitate those. 

The structure of the paper is as follows. We describe 
the setup in Section [TTl We discuss the Keldysh-action 
of the setup and introduce the counting field required 
for computing the statistics in Section IIIII We evaluate 
the FCS in Section IIVI and present the results for the 
photon Fano factor and big deviations from the equilib- 
rium. In Section |V] we give the interpretation in terms 
of bursts computing the partial rates of corresponding 
2fc-photon processes. We discuss two limiting cases of 
the FCS described in Section IVII We analyze the time- 
dependent fluctuations of radiation intensity in Section 
IVIII making use of the Keldysh propagator of the fields. 
In Section IVIIII we discuss the experimental significance 
of the frequency-resolved intensity correlations and quan- 
tify those. We conclude in Section HXl and give details of 
the field propagator in the Appendix. 



II. SETUP 

We concentrate on a setup similar to Ref. 0. In main, 
it comprises a Josephson junction biased by a d.c. volt- 
age source that is connected to a high-quality (that is, 
quality factor Q ^ 1) microwave resonator (Fig. [T]). We 
describe the resonator losses by the damping rate F. All 
photons leaving the resonator are absorbed by a (count- 
ing) detector. It is characterized by an efficiency /, a 
fraction of photons that are successfully counted. The 
impedance near the resonant frequency f2o reads 



Z{u! « Qq) 



-iv + T/2 



(1) 



V = uj — <^ VLq being the frequency mismatch. We will 
mostly concentrate on quantum limit of vanishing tem- 
perature UbT Ml. In this case, no photons come from 
the environment and the detector reading is the number 
of photons emitted from the resonator. The setup is char- 
acterized with a single quantum variable 4>{t) , related to 
the voltage across the inductor by means of Josephson 
relation = 2eV{t)/h. The superconducting phase dif- 
ference across the junction, is contributed by (j){t) 
and the voltage source, — (f) + 2eVbt. 

We will assume that the impedance far from the reso- 
nance, Zq, is sufficiently small at the quantum scale, that 
is, ZqGq ^ 1, Gq = e'^/Trh. Under this assumption, the 
junction is effectively in a low-impedance environment, 
and the contributions to the quantum fiuctuations of 
(j)(t) coming from frequencies far from VIq, d(j} yj ZqGq, 
can be safely neglected. Since practical impedances are 




FIG. 1: Setup. The Josephson junction with Josephson en- 
ergy Ej (cross in the Figure) is connected to a resonator repre- 
sented with an inductor and a capacitor. All resonator losses 
are absorbed by a detector and converted to a measurable 
signal with efSciency /. Right pane up: the impedance of 
the resonator in the vicinity of the resonant frequency Qq. 
Right pane down: the emission intensity from the resonator 
as function of Ej, Et corresponds to instability threshold. 



in the range of tens of Ohms, this assumption is well- 
justified. We stress that the assumption does not restrict 
the impedance near fioj — Zq{VLq/T) that can exceed 
the quantum scale at suflticiently big quality factors. This 
however is not needed for our approach to be valid: we 
only require Q ^ 1. 



III. KELDYSH ACTION 

Quantum dynamics of Josephson junction are well- 
explored-^^. The most general and adequate quantum 
description of the setup is provided^^ by a Keldysh- 
type path-integral over the variables (/>^(i), ± refers to 
the c-values of the quantum variable (t){t) at the for- 
ward(backward) part of the Keldysh contour. The "par- 
tition function" Z that is identically 1 in the traditional 
Keldysh approach is given by the path integral over the 
configurations of <i)^{t) weighted with the factor e''^, S 
being the quantum action expressed in Keldysh variables 
(from now on, we set h = 1) 



V[cj>+{t)]V[r{t)]e 



iS 



(2) 



The whole action is composed from Josephson and envi- 
ronmental part, S — Scnv + Sj. 

The action of the Josephson junction is simply given 
by its energy i?(0j) = —Ej cos{(j)j) and reads 



Sj 



Ej / dt [cos 0J (t) - COS (j)j (t)] 



(3) 



The superconducting phase difference across the junction 
is contributed by the d.c. bias voltage Vb applied, so we 
substitute 

<l>j{t)^^Vbt + <l>{t) (4) 
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into this part of the action. 

The action of the environment for a general frequency- 
dependent impedance Z{uj) reads 



with 



5*0: 



SttGc 



^ J2 (5) 



2n 



a.,l3=± 

with /dte*"V(0. and 



M{u) ^uj ^lmy(cj) 



(2n(cj) + 1) Rey(w) 



1 

-1 

1 -1 
-1 1 



Rey(w) 



-1 

1 



(6) 



where n{uj) = (exp(a;/fcBT) — gives the Bose- 

Einstein filhng factor at temperature T and the admit- 
tance Y{ll!) = Z^^{ljj). Variation of the action S with 
respect to 0+ — 0_ at 0+ « 0" f« (/) reproduces the 
"classical" equation of motion that disregards thermal 
and quantum fluctuations of cj), 



It: Ze 



2ei;jsin0(t) = (7) 



and is equivalent to condition of current conserva- 
tion. The admittance Y{uj) here determines the time- 
dependent response of current on voltage 0/2e. 

We specify to the case of a single resonance mode, such 
that the impedance near the resonant frequency ^Iq is 
given by Eq. [TJ To achieve the conditions of the paramet- 
ric resonance, we tune the d.c. bias voltage to Vb = hSl/e 
corresponding to the Josephson frequency 2VL close to 
the double of the resonant frequency fio- The detuning 
= ^ — VLq is assumed to be much smaller than fio- 
To implement this assumption, we introduce a slow com- 
plex variable <^(t), an amplitude of the resonant field, and 
express the original variable </> as 



(j,^{t) = 2nt + Re[e~*"V^(0] 



(8) 



thereby disregarding its Fourier components far from ±17. 
This is equivalent to a rotating-wave approximation. 

We substitute (f){t) to Eq. [3]in the form ([8]) and average 
it over the period of resonant oscillations to obtain local- 
in-timc action for the slow variable ^{1), 



Sj^jdt {Sj{ip+{t)) - Sj{ip-{t))) ; (9) 

Ej J2(|^|) . 2 
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(10) 



We also express the environment part of action in 
terms of the slow variable. 



(11) 



1 

-1 



nn + \ - 
-{nn + 1) no 



(12) 
(13) 



where we introduce integration over "low" frequencies v. 
The above expression can be rewritten in local-in-time 
form, that contains local time derivatives of the fields 
only. 



IGttGqZq 



dt (^ip~^* dt(p'^ — (p *dt(p 



/ +* + I -* - 



-r (no(^+*(^" + [nn + 1)(^"V^)) 



(14) 



This form of the (part of the ) action is proficient to es- 
tablish a relation with traditional optical techniques. If 
we rescale the variable 0, b = (f)/ {A^JGqZq), the rescaled 
variables 6*, h will provide path- integral representation of 
creation/annihilation operators h\h satisfying the stan- 
dard commutation relations. Since the action is local 
in time containing the derivatives only, the path-integral 
can be solved with an evolution equation. In case un- 
der consideration, this evolution equation is the Bloch 
equation in the rotating-wave approximation for density 
matrix in b\ b variables. It assumes the standard form 
implemented, for instance, in^. We do not outline this 
equation here since we proceed with a different method. 

Within this approximation, the " classical" equation ([7]) 
that corresponds to the saddle-point solution of the ac- 
tion can be written as 



dip 



<p + i{%nEjGQZo) 



p{{p*Y + v')- 



(15) 



a,l3=± 



(See e.gji^). The Eq. (fTS)) has stationary stable non- 
trivial solutions iy9 7^ 0, provided the Josephson energy 
exceeds a threshold Ej > Et, Et = fl/4nGQ\Z{n)\ = 
(r^ + 4:1^^Y^^/{8ttGqZo). These solutions give coherent 
emission at the frequency il. Below the threshold, quan- 
tum fluctuations enable emission of photon pairs result- 
ing in incoherent radiation with linewidth ~ F. 

We will restrict our consideration to the situation be- 
low the threshold. It is essential to note that the typical 
quantum fluctuation of if remains small below the thresh- 
old, {6(p{t))'^ ^ 1. This is guaranteed by the fact that 
the impedance Zq is small, {6ip{t))'^ ~ ZqGq. The fluc- 
tuations eventually grow at approaching the threshold. 
However, they will become of the order of 1 only in a 
close vicinity of the transition estimated as | Ej — i?t | ~ 
{ZoGQ)Et ^ Et. Therefore, almost everywhere below 
the threshold we may expand Sj in Taylor series in ip 
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keeping the leading quadratic term only, 

Sj = ^ jdt{{^^f-{v-f+c.c.) (16) 

We conclude that below the threshold the total action 
is quadratic describing Gaussian fluctuations of the field. 
One could get an impression of rather trivial statistics. 
Indeed, if we consider the statistics of the field itself, as 
it has been done in^'^ for a general linear electric circuit, 
we would end up with normal distributions. The point 
is that we are interested in the statistics of the photon 
flow, a variable that is quadratic in field. This leads to a 
non-trivial non-Gaussian statistics. 

Our goal is thus to describe the full counting statis- 
tics of photons emitted from the resonator. Most gen- 
eral characteristic function of these statistics is expressed 



Z{{xm ^ Tr 



Texp ( -I / dti{t)^ ) 



Texp -i / dti{t) 



(17) 



where T(T) denotes (anti)time ordering of the exponents, 
I = dtN is the operator of photon flow from the res- 
onator, N being the photon number opercitor, p—oo be- 
ing the density matrix . Indeed, expansion of (jl7p in 
powers of x(<) delivers the time-dependent correlators of 
the operators /. The characteristic function can be pre- 
sented by a path integral over the field configurations 
with the Keldysh action modified by the counting field 
x(i). (seei^ii^ for fermion casCjifl for photon case). With 
this, the only modified term in the action is the third one 
in , and the modification reads 



-(no + 1) nn + 



-{nn + 



l)e 



ix{t) 



-nne 



2. 



(18) 



This form of the modification is suggestive and can be 
derived heuristically. The fact that counting field enters 
the action in the form of exponents guarantees the inte- 
ger number of counts. If one rewrites the action in the 
form of master/Bloch equation for an extended density 
matri x ii ' ^ ' i ' " ' " ' ' ' !^ -, the modification concerns the terms that 
describe transitions with emission {T exp(ix){nn + 1)) or 
absorption (rexp(ix)nn) of a single photon, filling fac- 
tor of the environment photons entering the rates of these 
transitions in an expected way. 

The time-dependent counting field in the action is a 
parameter, that can be chosen at will. A common choice 
is a piecewise-constant x(^) = X with a time in- 

terval (0,r). Computed Z{x) becomes in this case the 
characteristic function of the probability distribution of 
emitting N photons within this time interval , 



ixN 



(19) 



and the cumulants of N are obtained via the differential 
relation 



((7V'"))=Cln(Z(x))lx=o. 



(20) 



In this work, we will concentrate on the low-frequency 
limit of the FCS assuming r to be much bigger than the 
typical waiting time of the phonon emission and disre- 
garding the contribution associated with the ends of the 
interval that does not depend on t. With this. 



(21) 



all information about the statistics being incorporated 
into a dimensionless function A(x)- The advantage of 
this assumption is that one can disregard the time- 
dependence of x{t) in the action that automates the eval- 
uation of the path integral. 

So far we have assumed an ideal efficiency of counting. 
If the statistics in this limit are known, one can easily 
obtain the results for any efficiency /. The method is to 
replace in all expressions for characteristic functions 



exp(ix) 1-^1 + /(exp(ix) - 1). 



(22) 



It is simple to justify this heuristically. One can split the 
whole damping rate F into undetectable losses Fi and 
losses detected, F2, / = F2/(Fi + F2). Fi and F2 both 
provide independent additive contributions to the action, 
and only the second one is modified with the counting 
field. 



IV. FULL COUNTING STATISTICS 

To represent the resulting action in a compact form, we 
introduce four independent scalar real fields correspond- 
ing to the complex fields ip for positive and negative v at 
forward/backward part of the contour. We group those 

in a 4- vector V'ly = (ft > j > (V-i/) > such that 

the modified action can be expressed compactly as a 4 x 4 
quadratic form in and tpl, 



dv 
2^ 



( M{v,x) A 

I A M^{-v,x) 



(23) 
(24) 



where 2x2 matrix M{v, x) is given by Eq. [T^ with the 
modification and 



A 



EV 



1 
1 



(25) 



where we have introduced a convenient dimensionless 
measure of Josephson energy E = 8ttGqZo{Ej /T). 

We take the path integral. Since it is Gaussian, the 
computation amounts to evaluation of the determinant 
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of the quadratic form. Since x is assumed to be time- 
independent, the quadratic form separates for each fre- 
quency. It is convenient to introduce discrete frequencies 
spaced with 2tt/t. Then, the integrals over the fields 
at each discrete frequency are Gaussian integrals propor- 
tional to the inverse of the determinant of matrix yl^. 
We transform the resulting product of determinants into 
the exponent of a sum and go to the continuous limit 
in this sum recovering the integral over the frequencies. 
The result in the integral form reads 



2{x) = cxp 



— In 

27r 



det jA^ix)) 
det(A,(x = 0)) 



(26) 



It is convenient to introduce dimensionless variables: 
i> — 1v jV ^ i>o ~ 2i/o/r, and a dimensionless parameter 
d = 1 + Vq — E^. The latter is important and enters most 
results presented below. The parameter d is positive, 
d = at the instability threshold, d = 1 + i>q > 1 at 
Ej = 0, this is, in the absence of the parametric driving. 
With this, the statistics are expressed in a simple integral 
form. 



A(X) 



^Infl + i^ 



Z(X) = E'nl (1 - + 
p{D) = D^ + i)^ 2{2-d) + d^ 



(27) 
(28) 



(29) 



We take the integral over the frequencies to arrive at 



A(X) 



-1 



This gives the FCS at arbitrary temperatures. The struc- 
ture of z{x) suggest that photons are emitted/absorbed 
in pairs (exp(±i2x) factors). Each emission/absorption 
probability is affected with the filling factors as ex- 
pected: absorption probability is proportional to Hq 
(two photons) , while emission is stimulated with a factor 
(1 -|- no)^). In the limit of small E <^ d one can expand 
A in terms of z to arrive at 

A(X) = -^^nie-'^'' - 1) -^{l + nnne''^ - 1) (31) 

This limit corresponds to the independent pair emis- 
sion/absorption acts, that can be thus regarded as un- 
correlated events. This gives Poissonian distribution of 
pair coimts. The rate of pair emission(absorption) is 
Te = T{E'^/4:d){l+nn)^){Ta = r{E'^/4:d)nl) and is much 
smaller than F under assumptions made. Upon increas- 
ing E, the correlations between the pair events set in. No 
event would take place at zero E. 

It may seem strange that at finite temperature no 
single-photon events are manifested in the FCS we 
present. Such events do take place, even in the absence of 
the parametric drive E: photons from the environment 
are randomly absorbed/emitted by/from the resonator. 



The point is that we concentrate here on the statis- 
tics in zero-frequency limit, and count all photons emit- 
ted/absorbed. In terms of cumulants of counts within a 
finite time interval r, we thus concentrate on the part 
of a cumulant that grows cx r. Such parts are absent 
for the single-photon statistics mentioned, and therefore 
these events do not contribute to the FCS we describe. 
An alternative way to understand this is to notice that 
without parametric drive the detector is in thermal equi- 
librium with the resonator. It is known that in this case 
it will not produce any (count) signal. 

From now on we will focus on the quantum limit, that 
is, on the case of vanishing temperature ksT <^ Ml, so 
that nn ^ 0. In this case, z{x) = E"^ {l - e^^x). This 
indicates that in this limit only pair emissions take place. 
The FCS expression reduces to 



A(X) 



+ +E^{l-e^^x) (32) 



This is one of the main results of this paper. In the 
limit of zero detuning d = 1 — E"^ and r — > cx3 this coin- 
cides with the results of Ref. j^. 



A. Average intensity and Fano factor 

Let us evaluate the first two moments of the statistics 
derived: the average intensity / = {N) /r and the in- 
tensity noise Si = {(N^)) /t. Expanding Eq. [32] in x 
gives 



(TV) 

T 



F dXix) 
2 ditx) 



(illustrated in Fig. [T]) and 



m) 



d'Xix) 



d{ixy 



AN) 



(4 + rf) 

d3 



(34) 



From this we can determine the Fano factor F — 
{(N^)) / (N) = Si /I. This number is significant in elec- 
tron or photon counting statistics giving an estimate of 
a number of particles that correlate with each other. We 
obtain 



F 



E' 



d2 



(35) 



We see that Fano factor is 2 in the limit of weak paramet- 
ric driving and diverges upon approaching the instability 
threshold d — >■ 0. It is instructive to re- write this ex- 
pression in terms of average number of photons in the 
resonator N = I/T — E"^ /2d and detuning vq. 



F ^2 + 2N [1 



A{2N + \) 



(36) 
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FIG. 2: Fano factor versus Josephson energy Ej. The en- 
ergy El corresponds to the instabihty threshold. Left paneh 
dependence on detuning uq. The curves from bottom to top 
correspond to = 0, 0.5, 1, 2, oo, sohd curves corresponding 
to the extreme values = Q,oo . Absolute efficiency / = 1 
is assumed. Right panel: Dependence of F on detection ef- 
ficiency f at uo — 0. The curves from top to bottom cor- 
respond to the efficiencies / = 1,0.5,0.25,0.1,0, solid curves 
corresponding to the extreme values / = 1, 0. 



The increased Fano factor is surely due to photon 
bunching. A naive picture of such bunching would pre- 
sume that N photons (in the resonator) stimulate emis- 
sion of one another, that is F on N. Hanbury-Brown- 
Twiss relation also supports such estimate. In the limit of 
large detuning ;>o — oo we indeed recover F = 2(1 -f N). 
However, generally it is not so: near the instability 
threshold {N ^ oo) F = 16^V(1 + i>^) ~ N^. The 
number of photons correlated exceeds by far the number 
of photons present in the resonator! This is specific for 
the parametric resonance. 

If only a fraction of the emitted photons is measured 
by the detector, the correlation decreases. With the aid 
of Eq. [22] the Fano factor measured can be expressed in 
terms of the Fano factor at absolute efficiency. 



F{f)^l + f{F-l) = {l + f)+E^ 



<P 



-f (37) 



It approaches 1 in the limit of small efficiency. 

The dependence of the Fano factor on the parametric 
drive is illustrated in Figure [21 



B. Big deviations 

The FCS expression (15^ can be employed to find with 
the exponential accuracy the probability of big deviations 
of / from its expectation value / (see e.g/^S,). To do this, 
one evaluates integral in Eq. [TH] ai N = It in the saddle- 
point approximation to obtain 



P{I) 



2tt ' 



-Mx) 



oc e 



-C(I) 



with 



£{x) = mm (jTj^f^ + 



(38) 



(39) 



T 1 r 



_l I S L 



//{r/2) 



FIG. 3: Probability of big deviations of J for £ = 0.7, uq = 0. 
(/ « 0.96(r/2)). Dashed line: Gaussian approximation valid 
for small deviations from the expectation value. 



For any FCS expression, £(/) = 0, and achieves a mini- 
mum there. The quadratic expansion near the minimum 
corresponds to a Gaussian distribution of small devia- 
tions from the expectation value. The probability of big 
deviations is not Gaussian although exponentially small. 

A typical dependence of In P on J is shown in Fig. [3| 
along with its Gaussian approximation. The probability 
is lower than the Gaussian approximation at / < / and 
higher otherwise. 

A feature worth discussing is that the probability to 
emit no photons (/ = 0) is finite and given by 



- — lnP = A(z(X)) = -l+\/l-f + V(f) (40) 

Rather counterintuitively, this probability remains finite 
even at approaching the threshold where / — > cxd. 



InP 




-1 



1 



1 



(41) 



Another feature worth discussing is the probability at 
I ^ I. The log of the probability appears to be propor- 
tional to /, 



— In P = /r/io 



(42) 



where — i/io gives the position of the singularity of X{x) 
in the plane of complex x- The singularity comes either 
from the inner or outer square root in ((321) ■ Owing to 
this, /xq exhibits a peculiarity (discontinuity of the second 
derivative) at d = 2 where the square roots merge into a 
1/4 singularity (Fig. [J) 



fJ-o 




(43) 
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FIG. 4: Probability of big deviations I I (EqgS} exhibits 
a peculiarity at d = 2. Left: The "transition" line d = 2 in 
the plane of parametric drive E and detuning vo separates 
the plane into the regions of small and large detuning. Right: 
The coefficient /io plotted along the dotted line in the left 
pane. To make the peculiarity visible, dashed curve gives the 
analytical continuation from d < 2 to d > 2. 



The condition d = 2 or, equivalently, E = ^v^ — \ gives 
thus a "transition" line that separates the parameter re- 
gions with large and zero detuning (Fig. 21). 

It is also possible to find the next-to-the-leading term 
in the asymptotic expression l^^ . an offset of linear 
asymptotics visible in Fig. |4l so the asymptotics become 



InP = /r/io - C; 



C 



7(1 - d/2) 



if d < 2 
if d > 2 



(44) 



Their dependence on E is illustrated in Fig. [SJ At small 
E, Ffe ~ E"^^ as expected for the rate of an event encom- 
passing 2fc photons. 

Note that the rates do not diverge at the threshold: 
rather, they saturate at finite value that decreases with 
increasing fc. To reconcile this with divergence of the 
radiation intensity at the threshold, let us determine the 
asymptotic behavior of Ffe in the limit of large k. The 
integral in Eq. 27] is contributed by minima of piy) , and 
can be approximated by a Gaussian integral. There is a 
single minimum at = if d < 2 and two minima at 

V 



±y/d — 2. The integration gives 



1 



fc3/2 



exp(-2/iofc)^o, 



(48) 



where 



d^+iE^ if d < 2 

if d>2 



(49) 



Comparing this with the probability of big deviations, we 
conclude that the big deviation is most likely a result of 
a single burst encompassing k = It photons during the 
observation interval. 

Near the threshold, these asymptotics read 



fc»i 



Et 



r exp {-k^ 
80F fc3/2 



(50) 



V. INTERPRETATION: BURSTS 

To understand better the FCS (l32|) . let us give an inter- 
pretation of these statistics. Let us note that the integral 
form of A(x), Eq. (1271) permits an expansion in powers 
of e"^^^, 

Jo 2^ [p{i>) + AE^) 
We rewrite it in the form 

p oo 

- 2 = E (exp(z2fcx) - 1) . (46) 
fe=i 

This suggests that photons are emitted in course of un- 
correlated events, bursts, each accompanying k photon 
pairs. The rate of a fc-burst is given by 

Analytical expressions for F^ become increasingly com- 
plicated with increasing k and we do not give them here. 



We see that at the threshold the relative probabilities 
of /c-bursts satisfy power law k~^/^. Although the prob- 
ability of big bursts is low, their contribution to the ra- 
diation intensity is high such that the average intensity 
diverges. Below the threshold, the power-law distribu- 
tion is cut at fc ~ TV^. This gives an estimate of the 
typical burst size contributing to the intensity, which is 
in agreement with an earlier estimation obtained from 
the Fano factor (Eq. 155)) . 

The fc-dependence of the rates at not-so-big k is illus- 
trated in Fig. [5j We see that the asymptotics is reached 
at rather low fc. 



VI. LIMITS 

We consider here two specific cases of the FCS under 
consideration: the limits of large detuning i7o ^ 1 (large 
d), and small d (the vicinity of the threshold). 



A. Large detuning 

In this case, one can assume c? ^ 1 almost everywhere 
in the sub-threshold region except a close vicinity {Et — 
E)/Et ~ of the threshold, and E < ^/d. In this case. 
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FIG. 5: Rates of the bursts. Right panel: The curves from 
top to bottom give the rates Fi through F5 at vq = ver- 
sus _E, i5 = 1 is the threshold . The rates are plotted as 
a function of Josephson energy scaled with respect to the 
energy Et, corresponding to the instability threshold. Left 
panel: The dependence of F^ (normalized on the power- 
law k~^''^) on the number of photon pairs in a burst. The 
curves from top to bottom correspond to Josephson energies 
Ej/Et = 1, 0.75, 0.5, 0.25, 0.1. Integer values of k are marked 
by squares. 



imation at time-scales exceeding /. We rearrange terms 
to arrive at 



A(X) 



-1 



+ ^l-f (l-yi-Szx Eyd^) (52) 

From this it is clear that small x ^ (P/4:E'^ eventually 
determine the statistics. Now we can expanding in d to 
recover a simpler FCS expression 



(53) 



To find the probability P{I) of big deviations, we take 
the integral in saddle point approximation 



P{N) - exp 



TTd 



(54) 



This form has been discussed previously in the context 
of photon counting statistics^*^. The expression for the 
probability is valid only if exponentially small, that is, 
for the observation times t ^ (Fd)"^). This suggests 
the relevance of a long time-scale ~ (Td)"-^ ^ F"^ in the 
vicinity of the threshold. 



the inner square root in Eq. [32] can be expanded in . 
The resulting FCS depends only on N and reads 

A(x) = -l + -^l + ^(e2'x-l) (51) 

It corresponds to the 'naive' estimation of the Fano factor 
F ~ 2{N -1-1). This form is very similar to FCS of inco- 
herent light with a Lorentz-shaped spectral intensit y^°i^^ 
with N replaced by the maximum filling factor of the 
photons in the light. The difference is that in our 
case the photons come in pairs rather than one-by-one 
(exp(ix) exp(i2x)). 

The origin of this similarity is understood if we con- 
sider the spectral intensity of the pairs emitted. This is 
given by inverse of p(i') and in the limit of large detuning 
consists of two narrow Lorentz-shaped lines centered at 
!> = ±^/d. In course of pair emission, each constitute 
of the pair appear in a separate line. Since the lines do 
not overlap, the photon bunching takes place separately 
within each line and has the same form as in the single- 
photon case. 

B. Vicinity of the threshold 

The distance to the threshold is parametrized by d <C 
1, d = precisely at the threshold. We need to expand 
A(x) in d. The formal expansion, however, does not work 
resulting in expressions that are singular at x = and 
therefore cannot be associated with any probability dis- 
tribution. To preserve analyticity in x, we need to explic- 
itly address small x ^ 1- To this end, we may expand 
exp(2«x) in X up to the first order. This disregards the 
discreteness of the photon flow, which is a valid approx- 



VII. TIME-DEPENDENT CORRELATIONS 

The burst interpretation outlined above would have 
been fine if the uncorrelated bursts could be regarded as 
instant events. In fact, the events are not instant: it takes 
time to emit k pairs composing a burst. If TV ~ 1 and 
/c ~ 1, this time is of the order of F^^. Close to threshold 
, iV 3> 1 and the typical waiting time between pair emis- 
sions is short, ~ ~ r~^. However, a typical burst 
in this case encompasses N'^ photons. This implies that 
the time required for a burst is actually long, ~ T~^N, 
in agreement with the final remark in Subsection IVIBI 
The bursts are thus overlapped in time. Moreover, even 
in the limit of N <^ 1 when the events are pair emis- 
sions that are well-separated in time, the constituents of 
the pair do not have to be emitted simultaneously. The 
low-frequency FCS computed does not provide direct in- 
formation about such time correlations. 

Here, we will investigate the time correlations restrict- 
ing to a simple case where we can proceed perturbatively. 
Let us choose the time-dependent counting field in the 
form 

X(i) = Xi0(^i, ^1 + dh) + X2e(i2, t2 + dt2), (55) 

where e{ta,tb) = Q{tb - t)e{t - ta). The counting field 
is thus piece-wise constant that is non-zero in two time 
intervals. If the duration of these time intervals is small, 
dti,dt2 <C the chance to have photon emissions 

within these intervals is small and can be computed per- 
turbatively. In this case, the unperturbed action corre- 
sponds to X = and we expand in terms of the pertur- 
bation 

^int^^ . / dt{exp{tx{t))-lW+*{t)^-{t) (56) 

16GoZn ./ 
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(see Eqs. [HI [T5] ) We assume that the time distance 
between the intervals r = t2 — ii ^ dti,2 is much bigger 
than the interval durations. 

Expansion of the cumulant- generating function \nZ{x) 
up to the second order gives 

InZ(x) =Ci(exp(ixi) - l)dti + C2(exp(ixi) - l)dti + 
Cii(exp(ixi)-l)'(dii)2 + 

C22(exp(iX2) - l)'(dt2)' + 

Ci2(r)(exp(ixi) - l)(exp(zx2) - 1) (57) 

It is clear that Ci^2 give a chance of photon emission in 
the intervals, so that Ci^2 = I- C12 is of interest for 
us since it gives correlations between the emissions sep- 
arated by time r: if an emission has occurred within 
the time interval {ti,ti + dti), this increases a chance 
of emission within (^2,^2 + ^^2)- We express this in- 
creased chance in terms of a time-dependent excess in- 
tensity /ox(t), Ci2(t) = //cx(t). 

From the other hand, the perturbations give 



C12 — 



210(GqZo 



:((0+(r)0-(r)*0+(O)0-(O)*)) (58) 



Since the fluctuating field is Gaussian, all correlators 
can be readily expressed in terms of the field propagator, 

G„^(i,t') = (CW#(i')), (59) 

that depends on time difference only, Gap{t,t') — 
Gaf3{t " t'). The quantity of interest is expressed as 



Ci 



((0+(t)r(O)*)(r(i)*0+(O))- 



210(GqZo)2 
(0+(t)0+(O))(r(i)*r(O)*)) 

r2 



'210(GqZo)2 

r2 



(G34(t)G2l(t)+G3l(0G24(0) 



16 



(A1+A2) 



(60) 



The evaluation of the propagator is straightforward 
but cumbersome, so we present the details in the Ap- 
pendix. 

We calculate the correlation function of photon emis- 
sion events separated by a time t. Two contributions to 
the correlator read {"f± 



d2 (1 - d) 

d2 (1 - d) 



7_e 



1± VT 

-7+r|T|/2 _ 



d) 



7+e 



-7-r|r|/2 



-7±r|T|/2 



^ r\r\ 
E^ 



The resulting excess intensity is therefore expressed as 



r E^d 

4jr~dj 



-7±r|t| 



2i>ge-r|«l 
E^d 



It is instructive to introduce the number of excess pho- 
tons nex(T) emitted within the time interval — |t|, |t| and 
obtained by the integration of the excess intensity lex (t) 
over the time, ncx(O) — 0, 



E^d 



2(1 -d) 



± 



7± 



(62) 



where the total number of excess photons ria 
7icx(oo) is related to the Fano factor (Eq. [55]) 

rioo = - 1. 



(63) 



In the limit of small parametric drive — >■ 0, rioo — 1- 
This implies that each photon correlates with strictly one 
extra photon forming a pair. The time-dependence of the 
correlations in this limit is given by 



riexir) = 1 - exp(-rr) 



(64) 



not depending on the detuning. 

In the vicinity of the threshold, the correlations are big 
and mainly build up at the slow time scale ~ (Td)^^, 

nex(T) = F (^1 - exp ) (65) 

The time-dependence of nc-x{T) is illustrated in Fig. 15] 



CD 




FIG. 6: Time correlations of emissions: number of excess 
photons emitted in the time interval (— r, r) provided the 
emission of a photon took place at t = 0. From upper- 
most to lowermost, the curves correspond to Josephson en- 
ergies E = 0.75,0.5,0.25,0.1 at zero detuning. The dotted 
line marks the value riox = 1, one excess photon. The curve 
at small E only exceeds 1 only slightly, manifesting the fact 
that emissions occur in pairs. Emissions of pair constituents 



(61) separated by time interval ~ F 
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VIII. FREQUENCY-RESOLVED 
CORRELATIONS 



Experiments on PCS in Josephson parametric ampli- 
fier are plausible but, as all experiments on PCS, are 
difficult and long, requiring long times of data accumu- 
lation and careful characterization of extrinsic noises in 
measurement setups. The first experiments would most 
likely concern intensity noise, the second cumulant of 
PCS. However even in this case the measurement may 
be difficult since the measured signal has to be amplified 
and the amplifier brings in a substantial extra noise. A 
common way to avoid such difficulties in the context of 
low-tcmperature measurement^^ is to split a noisy signal 
into two parts and amplify them by independent ampli- 
fiers. The cross-correlation of two outputs will not be 
affected by the amplifier noise. 

In our setup, it is convenient to split the signal in 
frequency domain. We introduce two detectors absorb- 
ing emitted photons with frequency-dependent efficien- 
cies f 1^2(1^) ■ This results in two intensity signals 



du 
2^ 



dl\ 



(66) 



dl / dv being intensity per frequency interval. Por our 
setup, the average intensity per frequency interval reads 
(see Eq. (EH) ) 



dl 



dv 4!>2 + (j>2 _ df ' 



(67) 



To describe the PCS of the two signals, we introduce two 
counting fields xi,2- With this, the x-dependent part of 
the action reads 



^=..^/^((e--l)A(.) + 



IQttGqZo J 27r 



(68) 



We need only the cross-correlation of the intensities gen- 
erally defined as 



JI2 = T 



(69) 



dxidx2 

at X1.2 ^ 0. Employing perturbations in xi,2 we find 

Si2^ l^^fi{^i)f{i^2)S{vuV2) (70) 

where the intensity correlator is expressed in terms of the 
field averages as 



28(7rGQZo) 



Expressing the correlator in terms of the field propagator, 
we find 



S{vi,V2) 



2T- 



[4^2 + (j>2 _ ^)2]2 

[E\l + E^ + vl - 2vo - vlf5{vi + V2) 

+iE^ 5{ui - V2)\ (72) 

This defines the general form of the spectral-resolved 
correlations below the instability threshold. The corre- 
lations are delta-functional and persist only for exactly 
equal or exactly opposite frequencies. This seems to nat- 
urally describe bunching of the photons in the same fre- 
quency mode as well as emission of pairs with frequen- 
cies opposite owing to energy conservation. However, 
delta-functional correlations are an artifact of Gaussian 
approximation: taking non-linearities into account would 
result in a smooth frequency dependence. Since we inte- 
grate over relatively wide frequency windows, the exact 
shape of the smoothed delta-functional peaks is not im- 
portant for us. 

Most comprehensive choice of the frequency-dependent 
efficiencies is as follows: 



/l = Q{v - UJs),f2 = Qi^s ~ v)- 



(73) 



The fist signal is thus collected from all frequencies above 
the separating frequency , while the second one comes 
from all frequencies below Wg. The dimensionless nor- 



(71) 



malized cross-correlation S12 = Syij y I\l2 is plotted in 
Pig. [7] versus E at zero detuning and for several values 
of Wg. At low E, the correlations are formed by emis- 
sion of photon pairs at opposite frequencies. At = 0, 
the numbers of photons emitted in two windows are pre- 
cisely the same, this results in ideal shot-noise correlation 
s\2 = 1. At Wg 7^ only a part of the pairs are separated 
into different windows, so the correlation is smaller. The 
cross-correlation grows with increasing E owing to pho- 
ton bunching. At Wg = 0, the cross-correlation diverges 
at the threshold. At Wg 7^ 0, the growth changes to de- 
crease and the normalized cross-correlation vanishes at 
the threshold. The reason for that is the narrowing of 
the spectral intensity upon approaching the threshold, 
so that the correlated emissions concentrate in one of the 
windows. 



IX. CONCLUSIONS 

To conclude, we have studied full counting statistics 
of Josephson junction circuit in the regime of parametric 
resonance. This is important in view of recent experi- 
ments that enable the detection of full power dissipated. 
We present the interpretation of statistics in terms of 
bursts of multiple-pairs of photons. We support this 
interpretation by investigating the time-dependent and 
frequency- resolved correlations. 

So far our results are restricted to the parameter re- 
gion below the threshold where the field correlations are 
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Sl2 




FIG. 7: Normalized cross-correlation of intensities in two fre- 
quency windows (Ea l73|) versus E [uq — 0). The separating 
frequency ujs takes values 0,0.1,0.2,0.3,0.4,0.5 F from the 
uppermost to lowermost curve. 



Gaussian. It is very interesting to address full counting 
statistics and time-dependent correlations in close vicin- 
ity of the instability threshold where the non-linear ef- 
fects are important. This will be the subject of future 
research. 
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Appendix: Propagator 

For the perturbative calculations presented in the main 
text, we need the propagator of the fields or, equiva- 
lently, fields ip, defined as 



G„^(t,t') = (CWV'/3(t')), 



(A.l) 



at X = 0. We rewrite the action at x = 0,no = with 
the aid of a dimensionless matrix A,, 



32nGQZo 4 7 27r 
/ a(i>,i>o) 



(A.2) 



iE 

-2 a(-z>,-£>o) -iE 

iE a(-z>,i>o) -2 

V ~iE a{i>,-Do) 

a{x,y) = 1 -i{x + y). 



The propagator in frequency domain is readily obtained 
by inverting A^. 



G{iy) = 2\GqZoT-'^A-^{v). 



The determinant of the action matrix 

det(i) = (l>2+72)(i>2+^2) 



(A.4) 



has four (generally complex) roots at dimensionless fre- 
quencies ±Z7±, where 7-1- = 1 ± vl — rf- 

The propagator in time-domain is obtained by the in- 
verse Fourier transform. We separate advanced {t > t') 
and retarded {t < t') part of the propagator. For ad- 
vanced part, 



Ga 



where 4x4 matrices G\ read: 



+ _ IGTrGgZo 



7+Vl - A 
with temporary notations 



Gn 


-£;2 


Gi3 


Gi4 


Gn 


-£;2 


Gi3 


Gi4 




•-^13 






^14 


•-^13 


-^2 





Gn = + 2 (1 - Wo) - d - iv^ 
G13 = — iE [ \/l — d — iva 



f Vi- 


1) 






/Gn 


E^ 


Gi3 


Gi4 


Gn 


E^ 


Gi3 


Gi4 




•-^13 


E^ 




\ GI4 


•-^13 


E^ 


^11 



Gr - 



7+^ 



with temporary notations 

G21 = + 2 (1 + if>o) fVl -d + ivQ 



(A.3) 



G31 = -iE{2~ iDo + Vl - rf) 
G41 = iE (yT~~d + Wo) 



(A.5) 



(A.' 



Gu = iE\2 + ivo 
and 

/ Gn £^2 ([j^g \ 
_ le^GqZo I Gn Gi3 Gi4 

7_Vl — « ^14 "-^13 ^ ^11 I 
\ Z72 / 

\ Ui4 Lti3 -C/ / 

with temporary notations 

Gn = - + 2 (1 - Wo) fVl + ii>o' 

Gi3 ^ - iE (VT- d + Wo , 

Gi4 = -iE{2 + Wo - Vl - 
For the retarded part, 

Gn = G+ e-^^(*'-*)r/2 ^- e-^"(*'-*)r/2 . (A., 
where 4x4 matrices G^ read: 



-^2 


-^2 


Lr4i 


Gil \ 


G21 


G21 


'-'SI 


G^i 


G31 


G31 


'-'21 


G^i 


G41 


G41 


-^2 





(A.9) 
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and 



jpz rpz 

IGttGqZq I G21 G21 G^i 



~ 7-^/^^ ^31 G31 g| g| I ' (^-^^^ 

\ G41 G41 S E 

with temporary notations 
G21 = -E^ + 2{1 + ivo) (vT^ - i!>o) 
G31 = iE{2- ivQ - \/r^) 
G41 = iiJ ^\/l — d — ivQ^ 
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